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A new approach to derivation of various effective-field approximations for lat- 
tice spin models in a unique way is presented. It is shown that it can give a 
number of methods, including the DMRG method, that can be used to find gener- 
ally inhomogeneous solutions of 2D classical lattice problems. A method, closely 
related to the DMRG method but without necessity to perform spin renormal- 
ization, is derived, yielding results practically not different from the DMRG ones. 
The matrix-product wave function of Rommer and Ostlund can be constructed 
from the output of the method. The computational costs of all the derived meth- 
ods are smaller than those of the DMRG. Most of the results are applicable to the 
ID quantum systems, as well. 



1 Introduction 

The rapid development of computational technique is accompanied by an in- 
creasing interest in numerical treatments of the problems of statistical mechanics. 
Besides Monte Carlo techniques, which simulate finite systems, the effective-field 
methods, which, in principle, treat infinite systems, became recently extremely 
popular. 

After traditional effective-field approaches, as molecular-field approximation 
and quasichemical (Bethe) approximation, the first sophisticated effective-field 
approach was Baxter's Corner-Matrix Method [1], which gave the possibility of 
systematic improvements of the approximation. This was also the feature of the 
cluster transfer-matrix method of the author [2], [7-9], which was used for an 
effective description of the incommensurate structures. Finally, Density Matrix 
Renormalization Group (DMRG) Method for calculation of the ground state of 
one-dimensional quantum models was developed by White [3, 4]. A little later it 
was applied also to the 2D classical lattice models by Nishino and its relation to 
the Baxter's method was pointed out [5, 6]. 

All these works are closely linked to Kramers- Wannier [10] suggestion of fac- 
torizing the wave function. 

The DMRG method started from the real-space renormalization approaches, 
where, calculating the trace of the Hamiltonian or summing up the partition func- 
tion, spin degrees of freedom where reduced by discarding eigenstates of the Hamil- 
tonian corresponding to the highest energy level (lowest eigenvalues of the transfer 
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matrix). It was shown by White that, instead of the Hamiltonian, the density ma- 
trix of the target state is more appropriate for this purpose. In this paper we 
show that if the summation of the partition function is done in a proper way no 
further spin reduction is necessary, and the renormalization process by discarding 
the eigenstates of the density matrix can be eliminated from the DMRG method 
at all. 

It is shown that DMRG method is, in fact, not based on repeated reduction, 
renormalization, and addition of new spin variables as in the real-space renor- 
malization group methods, but on a subsequent summation of spin variables in 
presence of effective fields when calculating the partition function (energy of the 
ground state in quantum case). 

The paper is organized as follows: In Section 2 a general method for derivation 
of various of effective-field approximation in a unique way is developed. It is used, 
in Section 3, for re-derivation of the cluster transfer-matrix method and DMRG 
method, and derivation a new effective-field method closely related to DMRG, 
but working more effectively. In section 4 the links between our and Rommer- 
Ostlund's [11, 12, 13] approach are shown, and the correct choice of the density 
matrix in the case of non-symmetric or non-Hermitian Hamiltonian of kinetic or 
quantum models is elucidated. 

2 Method 

The main task of statistical treatment of a classical system of spins <ii described 
by a Hamiltonian H(a{) is calculation of the partition function 

Z=J2^P[-pH(a i )}. 

For short-range interactions the Boltzmann's weight exp[— (3H{ai)] can be 
rewritten in the form of a product of functions % which depend on spin variables 
of few rows. The number of rows is given by the range of the interactions. 

exp[-/Wf(<7i)] = . . . 7J_i(Ai_i, . . . , A i _i +fc )7^(A i , . . . , A i+fc )7; + i(A i+ i, . . . , A i+ i +fc ) 

(1) 

where Aj = . . . cr^j-i, cr^j, cr^ j +1 . . . are row spin variables and k is the interaction 
range. Usually the functions T are written in a matrix form and called transfer 
matrices. We shall treat T as a function, but using the common phraseology we 
shall call it transfer matrix (T-matrix). 

For the case of nearest-neighbour interactions the elements of the matrix 
T((jj 5 i, . . . , (Ti >n ; cr i+11 , . . . , <7 i+ i jn ) are equal to the values of the function 
T(crj 5 i, . . . , a i>n , 0"i+i,i, . . . , <7j + i jn ) 
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If the summation of the Boltzmann's weight for a 2D system is performed row 
by row, the calculation of the partition function can be reformulated in a series of 
consecutive steps 

^ K+i • • • > K+k-i, A i+k )^ +k (A i+1 , . . . , A i+k ) = A i $^ A ,_ 1 (A i , . . . , A i+fc _i) 

$P(A i; . . . , A i+ fc_i)7^(Ai, Ai+i . . . , Ai+k-i, A i+ k) = \i$>f +1 (A i+ i, . . . , A i+ k) 

At 

(2) 

or the same in the presence of the functions $ from the opposite half-plane 

Ai+fe 

where i = 1, . . . , N, and the $'s at the left-hand sides of the equations are results 
of previous summations. 'U' and 'D' are for 'up' and 'down', respectively, to 
distinguish from the next step where the summation will be performed in the 
horizontal directions. 

Eqns. (2) and (3) are equivalent, but in the effective-field treatment they yield 
different results, as the functions $ provide effective fields only from below in the 
first case or from below and above in the second one. 

The functions $^ are eigenfunctions of the T-matrices only for the homogeneous 
solutions deeply in the bulk; for inhomogeneous ones all the functions are spatially 
dependent and the upper and lower functions $^ D ' ) are different from each other 
(also for symmetric T-matrices). They represent Boltzmann's weights of effective 
fields by which one half of the lattice acts onto the other one. Similarly, and 
^f%^ +k can be considered as Boltzmann's weights of one dimensional problems, 
and, assuming that the functions can be factorized, they may be again treated by 
the T-matrix technique; e.g. = [7 T i: j, where Tij = $>fj%j. 

Let us study a general one-dimensional problem with n row of spins. Then, the 
T-matrices and "eigenfunctions" \1/ are defined on finite clusters for finite-range 
interaction problems. Further, for the sake of simplicity, we shall consider only 
models with interactions within a square plaquette. 

E r (°Xi-l» • • -><7n,j-i;<ri,j, ■ • .,<Tn,j)*f (<?l,j, ■ ■ -,^n,j) = Aj ((7l j-i , . . .,<7 n ,j-l) 

Wij} 

E • • • ) 0"n,j-l)T , ((7l ) j_l, . . . ,CT nJ _i;CTi iJ -, . . .,<7„,j) = Xj^((Tl,j, ■ ■ -,<Jn,j) 

{<7i,j-l} 

(4) 



(3a) 
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Fig. 1. Decorated one-dimensional lattice. • - old spins aij, o - new spins nij. 

Having in mind that we have to perform summation in (2) and (3), we try, in 
an effective way, to reduce the number of rows, n, in the one-dimensional model 
by one. For that reason we double each column and insert between them new 
columns with n — 1 spins, as shown in Fig. 1. 

The Boltzmann's weight of the whole ID lattice instead of 

• • - T j-l T 3 T j + l ■ ■ ■ 

would be 

• • • Tj-i P j-i,L P j-i,R T jPj[L P j,R T j+i ■ ■ ■ 

where the functions P are the transfer matrices between the original columns of 
spins, whose number is doubled, and the new ones 

P jlh = pT (&l,j, • • • ) a n-l,j, Cnj'i 7Tl,j 5 • • • 7 ^n-l,j) 

P j,R = P (ni,j, • • • , TTn-lj; Cri,j, ■ ■ , <T'n-l,ji a 'n,j) 

Tj = T(a[j, . . . , v'n-ij, o'nj', + ■ • ■ , (Tn-l,j+l,&n,j+l) 

The newly introduced spins tv (o) divide the lattice into isolated blocks so 
that it can be easily summed up over the original spins a and a' (•). After the 
summation the number of new spins is one less per column and the number of spin 
degrees of freedom is reduced, i.e. the right-hand sides of Eqns. (2, 3) are found. 

Since all the statistical properties of the ID lattice can be determined from 
the functions $ J 7 they will remain unchanged after introduction of the new spin 
variables if the following equalities are satisfied: 
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P T T Pi R ^ R = ^ R ^ L P T t Pi R = # L 



(5) 



All our further considerations are based on this requirement. 



Clearly, the functions P should be constructed from the known values of the 
function \1>. There are many possibilities how to satisfy (5). Here, we discuss two 
choices of P: 



(0 



where 



. , <t„_i,j, a nJ ; mj, iTn-ij) = *f (aij, a n j)(p(aij, . . . , (T n -i,j)) 2 

(6) 



p(<Tl,j, ■ ■ .,<7 n -l,j) = ^ *i (°"1.J' • • ■> CT 'n-l,J> <T n,j)*?( <7 l,J> ■ • ■ ' a n-l,j , 0"n,j) 

This choice of P leads to the cluster transfer-matrix method [2]. 



Pjlhi^lJ: •••■> a l+l,ji • • • > a n,j'i 7Tl,j) • • • > Kn-k,j) ~ 

= ( Tl ^'' • ■ •>^.J'» (7 '+1.J' • • - a n,j)P~*(Tl,j,- ■ .,Tlj;TT 1+1 j . . . 7^21, j)' 

M 

n-k=2l 



where p 2 is the square root of the inverse to the matrix 
p(Ti,j,...,TLj;itL +1 j...7r 2 i,j) = 



(7) 



p is the density matrix corresponding to the functions tyj . In the formulas for 
Pr, in both cases, the indices R and L should be interchanged. 

The construction of both P's are based on the idea that \1/ in (5) meets another 
\1/ from the neighbouring P, forms the density function p or density matrix p, which 
are canceled by the inverse in P, and finally, it replaced by the same \& from the 
second P in (5). 
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The choice (ii) of the function P reduces the number of spins by k = n — 21. 
k can be equal to 1 only if n is odd. 

If k = 0, the choice (ii) gives the Density Matrix Renormalization Group 
(DMRG) method. Now, no reduction of spin variables takes place and it must 
be done artificially by discarding the spin degrees of freedom corresponding to 
one-half of the density matrix eigenvalues (the smallest ones). In this case the 
requirement (5) is not satisfied for low-order approximations. 

In the DMRG, one-half of eigenvectors of the density matrix, play the role of P. 
Their sign may be chosen arbitrarily, what sometimes leads to irregularities in the 
convergence of the iteration process. The choice (ii) this problems removes. (Note 
that P corresponds to the eigenvectors of the density matrix, which, in distinction 
to Eqn. 7, is summed up over first half of the spins in the argument of \I/'s.) 

3 Some methods derived from the requirement (5) 

Cluster transfer-matrix method 

The choice (i) of P can be successfully applied to the Eqns. (2) if the possibility 
of factorization of the functions ^f +k = [Jf $y +k j(a i+1 j, . . . , a i+k j; a i+1 j +1 , . . . , 
<j i+ kj+i) is assumed. From the fact that the interactions are of short range, 
the same factorization follows for the transfer matrix % = Y[j ■ ■ ■ , o~i+k,j', 

. . . , (Ji + k,j+i)- Then, the left-hand side of (2) can be expressed as a product 
Ylj ^Y+k = rij Inserting P^Pr (choice (i)) between every pair of T, and 
summing up over all original spins a, we obtain a function on n — 1 rows of spins 
7r, i.e. the function \&f +k _ 1 

$U = TT $ U 

j 

^2 • • • ' Ki+k-l,j, <ri+k,j)(p(ni,j, ■ ■, 7Ti+fc-l,j))~ 5 • 

°i + k,j + l 

■ 3>Y +k j(7r i+ ij, n i+ k-i,j, &i+k,j; n i+ ij+i, . . . , n i+ k-i,j+i, a i+ k,j+i)- 

' 1i,j{TTi,j, • • • , TTi+k-l,], <7i+k,j] ^ij+l, ■ ■ ■ , TTi+k-lj+l, <T i+ kJ + l)- 

■ *f +1 (7r i) j + i, . . . ,7ti +k - 1 j +1 ,a i+k j +1 ) ■ (p(n iJ+1 , . . . ,7r i+fc _i J+ i))~2 (8) 

where the functions \1/ are calculated from the Eqns. (4) . 

The accuracy of the method is given by the number of rows k that have to be 
as large as possible, despite of the fact that, for our short-range interactions, the 
function T and $ may, in principle, be defined only on a two and one row of spins, 
respectively. 
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The cluster transfer-matrix method can be used for calculation of spatially 
dependent properties of inhomogeneous solutions of lattice models. If the inhomo- 
geneities are controlled by the boundary conditions, the method 'sees' only three 
of them, e.g. the lower, right, and left boundaries, as it is based on one of Eqns. 
(2) and Eqns. (4). If the boundary conditions are given on sides of a rectangle, in- 
formation from one of the sides is missing, i.e. the cluster transfer matrix method 
can be used only to a somewhat limited class of inhomogeneous situations. For 
homogeneous solutions all the quantities in (8) are independent of their indices 
and after necessary number of iterations the bulk value of $ is obtained. 

The cluster transfer-matrix method is much less time consuming than the 
DMRG method as no matrix diagonalization is necessary. Nevertheless, the value 
of the critical temperature of the 2D Ising model is only slightly worse than the 
DMRG one for the same size of $ as shown in Table 1. 

On the other hand, if the cluster transfer matrix method is applied to Eq. (3), 
it can describe an arbitrary inhomogeneous situation, but the resulting critical 
temperature is for low-order approximations distinctly higher then the exact one 
and converges with increasing k to the exact value rather slowly. 

The method can be generalized to long-, but finite, range of interactions and 
modified in several ways. This is described in more details in [2], [7-9]. 

Modified DMRG method 

To derive the standard DMRG method, we have to factorize the function $ 
(for nearest-neighbour interactions defined on one row of spins), according to Bax- 
ter [1], by a matrix product of matrices 0cr J )( T j+1 (^j; indexed by the lattice 
spins to which they belong: = EL, fiaij^ij+^ij', The multi-value 
variables £ acquire m values, which control the order of the approximation. 

To provide a deeper insight in the further formulas we shall depict the resulting 
equations in a graphical way. 

Let us denote the matrix a (^j 1 ; Cj+i) graphically by The black 

dots denote the spins aj,aj + i and thick vertical lines the multi- value variables 



Then, $ can be expressed as a matrix product . . . { [[ 11 11 11 1. . .. 

A summation is expected over the neighbouring thick lines, which describe the 
same variable Similarly, the close-to-each-other dots denote the same spin a. 




■u 




Analogously, <E> is . . . 




Then, the left-hand side of Eq. (3a) 



can be depicted as 



7 



l + l 



w 



w 



w 



w 



w 



w 



w 



w 



(3b) 



where % is written as a product of plaquette Boltzmann's weights W: % = 



Inserting between each two T-matrices T 



W 



functions P7 and Pr of 



the type (ii) constructed from the T-matrix eigenf unctions \1/ R and \I/ L , reducing 
the number of the rows and columns in the density matrix by one-half, summing 
up over all the spin variables a and £ in the row i, the standard DMRG method 
is obtained. Nevertheless, our approach, comparing with the standard one, shows 
more clearly how to to construct the density matrix in the case of non-homogeneous 
solution or non-symmetric T-matrix. 

In the one-dimensional system (3) even number of spins in each column is used: 
2 original spins of the model used in the Boltzmann's weight W and two m-valued 
spins £. To get odd number of spins, we must omit one of the original spins, 
i.e. also the weight W. This step would seemingly lead to omission of the whole 
information about the Hamiltonian of the system. To avoid this, we discard all 
W's but one. 



3 3+1 



i + 1 



W 



(9) 



In the diagram there appeared two new symbols , , which, as it will be seen, 
denote P^(i + l,j) and P^(i,j + 1). 

Summing over sites k < j in the row i, we get ^ .(£?., a^j, £R-)> and over the 
site j + 1 in the row i and all sites in the row i + we get , a^j , ) . From 

this two functions the density matrix is constructed according to the rule (ii) , and 
the functions P^(i,j) and Pn(z,j) are calculated. The density matrix is generally 
non-symmetric. For the homogeneous phases \E fR and \& L are the eigenvectors of 



the transfer matrix T = 



they are equal to each other, and the density 
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matrix is symmetric. The number of spins in \I/'s is three, and, after inserting 
-Fr(*)J) an d P^(i,j) between the matrices at the site their number will be 

reduced. 

The function P^(i,j) should be remembered for the calculation at the site 
— l and Pr(i, j), together with Pjj + 1) from the previous step, is used for 
calculation of the new value of 



3 



1 



n 



(10a) 



or 



■ W((Tij,(Tij + i, CTi+l,j, cr i+1 j +1 ) 



,3 

E 



(10b) 



cU cD eU tD 



The result of this summation is 



i + 1 



ui 



UUUIXl 



uuijnnn 
nnni 



(ii) 



In the next step a new Boltzmann's weight W can be created from dP- , • • 
between the columns j — 1 and j by a procedure inverse to (10), and the sweep to 
the left can be continued. The sweep to the right is performed analogously. 

When the steady state (fixed point) is reached, the sweeps leave the functions 
\&l and \1/r unchanged at any site of the lattice, and the calculations, in the 
framework of the approximation given by the number m of values of the variable £ 
in the matrices 0, are exact. The Boltzmann's weights can be created between two 
arbitrary rows so that the sweeps can be performed in the horizontal and vertical 
directions. Finally, after sufficient number of sweeps, the values of Boltzmann's 
weights of effective fields -Pl(r) an d U( ^ D ' ) are obtained at all the sites of the finite 
lattice, from which all interesting average quantities can be calculated. 

The calculation for homogeneous case is very simple. It consists of three steps 
iteratively performed until the fixed point is reached. 
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1. The ground state eigenfunction \1/ of the transfer matrix T 



is 



calculated. The functions $ j || J are assumed to be known from the previous 
step. 

2. Knowing the function P is found from (7). 

3. The new value of $ for the step 1 is calculated from Eqn. (10). 

As for the homogeneous case it is no need to distinguish between right, left 
upper, and lower functions, the indices R, L, U, D are ommited. 

In the effective-field methods spontaneous symmetry-breaking occurs, as a rule. 
Thus, the critical temperature is given by the onset of non-zero magnetization. 
Actually, no calculation of magnetization is needed, e.g. the difference 5>n(£, £') — 

(£'£') a ^ so snows if the system is in ferromagnetic or paramagnetic state. 

The described method is much faster than the standard DMRG method, as 
for the same accuracy the functions ^> acquire only one-half of the values of $ 
in DMRG and the density matrix has one-quarter of elements. This statement 
applies for the improved DMRG algorithms [14], as well, as it can be used for 
calculation of the superblock (T) wave function also in our case. 

However, far faster than the DMRG is the cluster transfer-matrix method, 
where instead of density matrix p a density function p is used, and the matrix 
inversion or matrix diagonalization is replaced by division. 

For comparison, the values of critical interaction constants (inverse of the tem- 
peratures) for all the above-mentioned methods in the case of the homogeneous 
2D Ising model are given in Table 1. 



Table 1. Critical values of the interaction constant K c of the Ising model on the 
2D square lattice. 



Method 


Number of spin-degrees 
of freedom in \1/ 


K c 


Cluster matrix method applied to Eq. (3) 

(fastest) 


2 10 = 1024 
(no diagonalization) 


0.4347 


Cluster matrix method applied to Eq. (2) 

(fastest) 


2 10 =1024 
(no diagonalization) 


0.44111 


Our modification of DMRG 


to = 16 2m 2 = 512 


0.44050 


Our modification of DMRG 


to = 22 2to 2 = 968 


0.44052 


Standard DMRG 
(slowest) 


to = 16 4m 2 = 1024 


0.44050 


Exact value 




0.44068 



Since the methods are based on an approximate summation of the partition 
function all interesting thermally averaged physical quantities can be found as 



10 



functions of coordinates. It was found that the simplified DMRG without renor- 
malization yields practically the same values as standard DMRG also for other 
relevant physical quantities. 

4 Discussion and conclusion 

For re-derivation of DMRG and derivation of the modified DMRG, the Bax- 
ter approximate expression in the form of matrix product for the wave function 
was used. Similar factorization was applied by Rommer and Ostlund (RO) [11] 
for derivation of their variational approach. Nevertheless, the matrices in both 
approaches are different. Baxter's matrices are indexed by two site spins while 
Rommer and Ostlund's only by one spin. They represent the functions P in our 
approach. In formula (9) they appear as vertical lines with one dot. It is easily seen 
that if the step in (9) is high of several rows the wave function \& in the vertical di- 
rection would be a product of RO matrices. The step from (9) to (11) may be taken 
not only as a calculation of (j)f +1 j. i+1 j +1 from but also as of P^(i,j) from 

P^(i,j + 1) or iterative calculation of RO eigenfunctions in horizontal directions. 
This procedure would lead to the minimum of the free energy. As formation of 
arbitrary steps is, in our approach, exact at fixed point, we see that we have to 
obtain the same results as RO method, if treated exactly for arbitrary m as it 
was done here. Our method and Rommer-Ostlund variational approach are com- 
plementary, but their results slightly differs from the standard DMRG ones, what 
can be seen, e.g. from the fact that XIo-^d -Pr(£ D > a i £ U )Pl(£ D i a , ?' U ) = Sgu^v is 
satisfied only for two first methods and not for DMRG. 

A direct application of RO factorization in the general approach described in 
Section 2 would lead to necessity of singular decomposition of matrices, what we 
want to avoid in this paper. We established a connection between Baxter's and 
RO factorization. The possibility of wave-function factorization in Baxter's way 
was assumed and RO factorization derived from it. 

We claim that our method avoids the renormalization of the spin variables. 
Here, the notion of renormalization is taken in the sense of real-space quantum 
renormalization group, where the states corresponding to some eigenvalues of cer- 
tain matrix are discarded and the spin variables are mixed together and lose their 
identity. This is not done in this paper. In both choices of P (i,ii) many 5-functions 
appear so that most of the new and old spins are really identical. However, the new 
spins 7T/+1J . . . 7T2ij in (7) are different from the old ones and may be considered 
in another sense as renormalized. 

The requirement (5) applies also to one-dimensional quantum and kinetic prob- 
lems. The matrices P^ and Pr have to be inserted between \1> and H in the 
expression ^ + H^, where H is the Hamiltonian of the superblock in the White's 
approach [2]. Now, from (5) it is clear that, for non-hermitian H, the left and 
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right reduction matrices (here P's) should be both equal and non-symmetric, or 
different and symmetric, constructed from the left and right eigenvectors, respec- 
tively. The fact that the standard derivation of DMRG does not give transparent 
prescription for constructing the density matrix in non-symmetric case is seen in 
[15], where instead of the correct choice [16, 17] of the non-symmetric density 
matrix is not used the second possibility of right and left density matrix at the 
left- and right-hand side of the superblock operator, respectively, suggested by us, 
but the incorrect average of the left and right density matrices, or only the right 
density matrix, or the density matrix with mixed terms. 

This approach enables effective numerical treatment of one-dimensional non- 
equilibrium kinetic models and quantum systems approximatively solved by ana- 
lytical methods [18]. 

It was proposed a new method for derivation of different types of effective-field 
approximate methods for calculation of thermal averages of physical quantities 
and thermodynamic functions of inhomogeneous classical 2D spin lattices and 
zero-temperature properties of ID quantum lattice system. 

The effective fields appear after summation over all spins except a square pla- 
quette. In (10) they are applied from below, left, and right and are represented by 
(jp , , and Pr, respectively. The effective field from above affects the calculation 
indirectly, through and consequently P. The effective fields should be really 
perceived in a generalized sense - while in the cluster transfer matrix method their 
Boltzmann's weights are positive in the DMRG method, they may be negative. 
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